Sex & statistical power - simulations

Code

library(ggplot2)
library(pander)
library(emmeans)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ lubridate 1.9.2     ✔ tibble    3.1.8
✔ purrr     1.0.1     ✔ tidyr     1.3.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

library(here)

here() starts at /Users/szymek/Library/CloudStorage/Dropbox/#SCIENCE/R/230713_power_sex_sim

Code

library(ggpubr)
library(ggridges)

Analysis preps (colours, options, auxiliary functions)

Code

cbb_palette <- c(
  "#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442",
  "#0072B2", "#D55E00", "#CC79A7"
)

## Avoid table wrapping
options(width = 800)
set.seed(999)

source(here("PhillipsKarpEtal_data/funs.r"))

Below results repeat simulations from the paper of Phillips et al. (2018) PLoS Biol 21(6): e3002129. Original results simplified by removing the post-hoc and the (pooled) t-test methods of analysing data. Generation of simulated data was modified to use a formal linear model defined as:

\[\mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{e}\]

with \(\mathbf{b}=(\mu,\beta_{Treat},\beta_{Sex},\beta_{Treat:Sex})\) and \(\mathbf{e} \sim \mathcal{N}(\mathbf{0}, \mathbf{s})\), with \(\mathbf{s}=\sigma\sqrt{\mathbf{1}+\mathbf{X}_{Sex}\cdot\alpha}\) (\(\alpha\) quantifies the amount of heteroscedasticity in sex variances, with SD \(\sigma_{m}=\sigma_{f}\sqrt{1+\alpha}\), or male variance being \(\alpha\cdot100\%\) greater than female). Regression coefficients (\(\mathbf{b}\)) are compatible with contrasts (sex/treatment effects) used in Phillips et al. (2023).

1 Scenarios from Fig. 3

1.1 From original paper

Assuming homoscedasticity of sex variances

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0

all_effs_sex_none <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_none <- all_effs_sex_none %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

calc_power_sex_none %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1

all_effs_sex_small <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_small <- all_effs_sex_small %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none"))

calc_power_sex_small %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1

all_effs_sex_large <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large <- all_effs_sex_large %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("none"))


calc_power_sex_large %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Plots combining all scenarios from Fig. 1

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_none, calc_power_sex_small,
                               calc_power_sex_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(sex_eff_size = factor(sex_eff_size,
                               levels = c("none", "small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = sex_eff_size,
             group = interaction(sex_eff_size, method))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = method)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Sex effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Stat method")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  # filter(!Effect %in% c("T-Test Treatment Eff")) %>%
  mutate(sex_eff_size = factor(sex_eff_size,
                               levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
             group = interaction(sex_eff_size, Effect))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sex_eff_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sex effect size") + ylim(0, 1)

1.2 Heteroscedastic sex variances

We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.

Low heteroscedasticity

Code

alpha <- 1

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_none_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_none_lowh <- all_effs_sex_none_lowh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_none_lowh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_small_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_small_lowh <- all_effs_sex_small_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
  mutate(heterosc = rep("low"))

calc_power_sex_small_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity

Code

alpha <- 2

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_none_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_none_largeh <- all_effs_sex_none_largeh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_none_largeh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_small_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_small_largeh <- all_effs_sex_small_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_small_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

It should be of note that evolutionary biology has seen heteroscedasticities exceeding 100% in variance ration between sexes. Combined results below show the impact of heteroscedasticity.

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_small,
                               calc_power_sex_large,
                               calc_power_sex_small_lowh,
                               calc_power_sex_large_lowh,
                               calc_power_sex_small_largeh,
                               calc_power_sex_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(sex_eff_size = factor(sex_eff_size,
                               levels = c("small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = sex_eff_size,
             group = interaction(sex_eff_size, heterosc))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Sex effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity")  + ylim(0, 1)

Code

plotA <- all_sex_effs_sex_only %>%
  filter(sex_eff_size %in% c("large")) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)

plotA

2 Scenarios from Fig. 4 - interaction of magnitude

2.1 From original paper

Below simulations were only slightly modified. We assumed a moderate (+0.5) overall effect of sex and a varying interaction increasing the magnitude of sex difference in the treated group (+0, +0.5, +1).

Assuming homoscedasticity of sex variances

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0

all_effs_sex_ix_none <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1

all_effs_sex_ix_small <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1

all_effs_sex_ix_large <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("large"))


calc_power_sex_ix_large %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Plots combining all scenarios from Fig. 1

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_none, calc_power_sex_ix_small,
                               calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, method))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = method)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Stat method")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  # filter(!Effect %in% c("T-Test Treatment Eff")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
             group = interaction(ix_eff_size, Effect))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = ix_eff_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Ix effect size") + ylim(0, 1)

2.2 Heteroscedastic sex variances

We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.

Low heteroscedasticity

Code

alpha <- 1

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_lowh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_ix_small_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
  mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1+alpha

all_effs_sex_ix_large_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity

Code

alpha <- 2

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_largeh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_ix_small_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1+alpha

all_effs_sex_ix_large_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Power consequences of heteroscedasticities with varying (magnitude-only) interaction.

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_small,
                               calc_power_sex_ix_large,
                               calc_power_sex_ix_small_lowh,
                               calc_power_sex_ix_large_lowh,
                               calc_power_sex_ix_small_largeh,
                               calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, heterosc))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity")  + ylim(0, 1)

Code

plotB <- all_sex_effs_sex_only %>%
  filter(ix_eff_size %in% c("large")) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)

plotB

3 Scenarios from Fig. 5 - interaction of direction

3.1 From original paper

Below simulations were only slightly modified. We assumed a strong (-1) overall effect of sex and a varying interaction resulting of crossing of sex effects on the treatment effect (0, -1, -2).

Assuming homoscedasticity of sex variances

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0

all_effs_sex_ix_none <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1

all_effs_sex_ix_small <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1

all_effs_sex_ix_large <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("large"))


calc_power_sex_ix_large %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Plots combining all scenarios from Fig. 1

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_none, calc_power_sex_ix_small,
                               calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, method))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = method)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Stat method")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  # filter(!Effect %in% c("T-Test Treatment Eff")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
             group = interaction(ix_eff_size, Effect))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = ix_eff_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Ix effect size") + ylim(0, 1)

3.2 Heteroscedastic sex variances

We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.

Low heteroscedasticity

Code

alpha <- 1

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_lowh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1+alpha

all_effs_sex_ix_small_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
  mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity

Code

alpha <- 2

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_largeh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1+alpha

all_effs_sex_ix_small_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Power consequences of heteroscedasticities with varying (direction-only) interaction. Power anomalies indicate pretty low inferential abilities of such models at low sample sizes.

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_small,
                               calc_power_sex_ix_large,
                               calc_power_sex_ix_small_lowh,
                               calc_power_sex_ix_large_lowh,
                               calc_power_sex_ix_small_largeh,
                               calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, heterosc))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  filter(ix_eff_size %in% c("large")) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)

4 Impact of sample size on power

4.1 Large sex eff, no ix

Below simulation code generates power estimates for the following sample size scenarios: 5M+5F in each treatment group (original scenario, totql sample 20 individuals), 10M+10F (N = 40), 50M+50F (N = 200).

First - low heteroscedasticity:

5M+5F

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_lowh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh_5 <- all_effs_sex_large_lowh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
  mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_large_lowh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh_10 <- all_effs_sex_large_lowh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_large_lowh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh_50 <- all_effs_sex_large_lowh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity:

5M+5F

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_largeh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh_5 <- all_effs_sex_large_largeh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_large_largeh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh_10 <- all_effs_sex_large_largeh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_large_largeh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh_50 <- all_effs_sex_large_largeh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

4.2 Large sex effect, crossed reaction norms (interaction)

Identical sample size scenarios.

First - low heteroscedasticity:

5M+5F

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_lowh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh_5 <- all_effs_sex_ix_large_lowh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_ix_large_lowh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_ix_large_lowh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh_10 <- all_effs_sex_ix_large_lowh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_ix_large_lowh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_ix_large_lowh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh_50 <- all_effs_sex_ix_large_lowh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_ix_large_lowh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity:

5M+5F

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_largeh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh_5 <- all_effs_sex_ix_large_largeh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_ix_large_largeh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_ix_large_largeh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh_10 <- all_effs_sex_ix_large_largeh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_ix_large_largeh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_ix_large_largeh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh_50 <- all_effs_sex_ix_large_largeh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_ix_large_largeh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

4.3 Summaries for sample size scenarios

No interaction, low (left) and large (right) heteroscedasticity:

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_lowh_5,
    calc_power_sex_large_lowh_10,
    calc_power_sex_large_lowh_50,
    calc_power_sex_large_largeh_5,
    calc_power_sex_large_largeh_10,
    calc_power_sex_large_largeh_50
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot1 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("low")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)


plot2 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("large")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)

ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

Interaction of effect direction present, low (left) and high (right) heteroscedasticity:

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_ix_large_lowh_5,
    calc_power_sex_ix_large_lowh_10,
    calc_power_sex_ix_large_lowh_50,
    calc_power_sex_ix_large_largeh_5,
    calc_power_sex_ix_large_largeh_10,
    calc_power_sex_ix_large_largeh_50
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot1 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("low")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)


plot2 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("large")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)

ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

5 Simulating sample size

5.1 Simulations assuming target power = 80%

Scenario 1: strong sex effect, none to large heteroscedasticity. Assumed power 80% for treatment effect.

Code

alpha <- 0

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_noh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_noh,
    calc_power_sex_large_lowh,
    calc_power_sex_large_largeh
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) + coord_trans(y = "log10") +
  theme(text = element_text(size=16))

plot1

Code

plot1.1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
  xlim(0.3, 0.7) + ylim(5, 70)

plot1.1

Warning: Removed 36 rows containing missing values (`geom_line()`).

Code

alpha <- 0

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_large_noh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_noh,
    calc_power_sex_large_lowh,
    calc_power_sex_large_largeh
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot2 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) + coord_trans(y = "log10") +
  theme(text = element_text(size=16))

plot2

Code

plot2.1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
  xlim(0.3, 0.7) + ylim(5, 25)

plot2.1

Warning: Removed 36 rows containing missing values (`geom_line()`).

Code

alpha <- 0

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha

all_effs_sex_large_noh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_noh,
    calc_power_sex_large_lowh,
    calc_power_sex_large_largeh
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot3 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) + coord_trans(y = "log10") +
  theme(text = element_text(size=16))

plot3

Code

plot3.1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
  xlim(0.3, 0.7) + ylim(5, 2200)

plot3.1

Warning: Removed 36 rows containing missing values (`geom_line()`).

5.2 Combined plots

(In all below examples the y axis (sample size) indicates the number of individuals per group and sex - hence the total sample size is 4 times the number indicated on the y axis.)

Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for full range of effect sizes and on log-scale for sample size:

Code

comb1 <- ggarrange(plot1, plot2, plot3, ncol = 1, nrow = 3, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C"))

comb1

Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for moderate-sized effect sizes using identity scale for sample size:

Code

comb2 <- ggarrange(plot1.1, plot2.1, plot3.1, ncol = 1, nrow = 3, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C"))

Warning: Removed 36 rows containing missing values (`geom_line()`).
Removed 36 rows containing missing values (`geom_line()`).
Removed 36 rows containing missing values (`geom_line()`).
Removed 36 rows containing missing values (`geom_line()`).

Code

comb2

Combined plots demonstrating the impact of heteroscedasticity on power for a no-interaction case (A) and interaction-of-magnitude case (B; stronger effect in the more variable sex):

Code

comb3 <- ggarrange(plotA, plotB, ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom", labels = c("A", "B"))

comb3

6 Source of heterogeneity

Below plots show hypothetical scenarios when heterogeneity is a natural outcome of the underlying statistical process. Diagonal lines connect averages belonging to one treatment group (control vs. experimental). Thick vertical lines show SD of each subgroup.

6.1 Poisson distribution

Simulation draws from \(Pois(\lambda_i)\) assumes sex/treatment specific averages (\(\lambda_i\)) - see lambda1...lambda4 below.

Code

N <- 500
lambda1 <- 3
lambda2 <- 9
lambda3 <- 8
lambda4 <- 15
toy_poisson <- data.frame(
    Treatment = c(rep("Ctrl", 2*N), rep("Exp", 2*N)),
    Sex = c(rep("F", N), rep("M", N), rep("F", N), rep("M", N)),
    Response = c(rpois(N, lambda1), rpois(N, lambda2), rpois(N, lambda3), rpois(N, lambda4))
)
toy_poisson$Group <- interaction(toy_poisson$Sex, toy_poisson$Treatment)

Code

plot1 <- ggplot(data = toy_poisson, aes(x = Response, y = Group)) +
    geom_density_ridges(aes(fill = Sex, scale = 0.8),
        stat = "binline", bins = 15, alpha = 0.6, draw_baseline = FALSE
    ) +
    theme_ridges() +
    geom_point(size = 4, stat = "summary", fun = "mean") +
    stat_summary(
        fun = "mean", geom = "line", linewidth = 1,
        aes(group = Treatment, linetype = Treatment)
    ) +
    stat_summary(
        fun.max = function(x) mean(x) + sd(x),
        fun.min = function(x) mean(x) - sd(x), geom = "errorbar",
        linewidth = 1.5, width = 0
    ) +
    theme_bw() +
    coord_flip() +
    theme(text = element_text(size = 16))

plot1

6.2 Log-normal distribution

Simulation samples from \(Lognormal(\mu_i, \sigma^2=0.1)\) with sex/treatment specific means (mu1...mu4) below.

Code

N <- 500
mu1 <- 0.01
mu2 <- 0.5
mu3 <- 0.4
mu4 <- 0.8
sd <- 0.31
toy_lognorm <- data.frame(
    Treatment = c(rep("Ctrl", 2*N), rep("Exp", 2*N)),
    Sex = c(rep("F", N), rep("M", N), rep("F", N), rep("M", N)),
    Response = c(rlnorm(N, mu1, sd), rlnorm(N, mu2, sd), rlnorm(N, mu3, sd), rlnorm(N, mu4, sd))
)
toy_lognorm$Group <- interaction(toy_poisson$Sex, toy_poisson$Treatment)

Code

plot2 <- ggplot(data = toy_lognorm, aes(x = Response, y = Group)) +
    geom_density_ridges(aes(fill = Sex, scale = 0.8),
        stat = "binline", bins = 15, alpha = 0.6, draw_baseline = FALSE
    ) +
    theme_ridges() +
    geom_point(size = 4, stat = "summary", fun = "mean") +
    stat_summary(
        fun = "mean", geom = "line", linewidth = 1,
        aes(group = Treatment, linetype = Treatment)
    ) +
    stat_summary(
        fun.max = function(x) mean(x) + sd(x),
        fun.min = function(x) mean(x) - sd(x), geom = "errorbar",
        linewidth = 1.5, width = 0
    ) +
    theme_bw() +
    coord_flip() +
    theme(text = element_text(size = 16))
plot2

6.3 Combined plot

Code

#| fig-label: fig:source_hetero1
#| fig-cap: Simulated sex-by-treatment data using a Poisson (A) and log-normal (B) distribution. Diagonal lines connect averages belonging to one treatment group (control vs. experimental). Thick vertical lines show SD of each subgroup.

ggarrange(plot1, plot2, ncol = 2, nrow = 1,
common.legend = TRUE, legend = "bottom", labels = c("A", "B"))

6.4 Non-linear response in sexes

Heterogeneity in variances may also arise if the underlying phenotypic mapping function is non linear and has saturating/desaturating character in a range modified by sex dimorphism and/or experimental treatment.

Here we assume a logistic function:

\[f(x)=\frac{L}{e^{-k(x-x_0)}}\]

with parameter values: L = 3.5, k = 1, x0 = 1.

Code

mapping_logistic <- function(x, L = 3.5, x0 = 1, k = 1) L/(1+exp(-1*k*(x - x0)))

Using this function, we can generate data assuming some treatment and sex differences:

Code

SD_N <- 0.5
mu_env1 <- 1
mu_env2 <- 2
mu_env3 <- 3
mu_env4 <- 4.5
toy_nonlin <- data.frame(
    Treatment = c(rep("Ctrl", 2 * N), rep("Exp", 2 * N)),
    Sex = c(rep("F", N), rep("M", N), rep("F", N), rep("M", N)),
    Env = c(
        rnorm(N, mu_env1, SD_N), rnorm(N, mu_env2, SD_N),
        rnorm(N, mu_env3, SD_N), rnorm(N, mu_env4, SD_N)
    )
)
toy_nonlin$Response <- mapping_logistic(toy_nonlin$Env)

Code

#| fig-label: source_hetero2
#| fig-cap: Heteroscedasticity between sexes (colours) and treatment groups (solid/dashed lines) resulting from a non-linear (here - logistic) mapping of some manipulated environemntal variable into the measured response. Thin dotted line - underlying mapping function. Points - realisations of mapping for simulated environmental data. Histograms relate to observed response measurements.


ggplot() +
    geom_point(
        data = toy_nonlin,
        aes(y = Env, x = Response, colour = Sex),
        pch = 1, size = 3, alpha = 0.8
    ) +
    geom_line(
        colour = "black", linetype = 2,
        aes(y = seq(-1, 8, 0.1), x = mapping_logistic(seq(-1, 8, 0.1)))
    ) +
    geom_density_ridges(
        data = toy_nonlin %>%
            group_by(Sex, Treatment) %>%
            mutate(meanEnv = mean(Env)),
        aes(x = Response, group = meanEnv, y = meanEnv, fill = Sex),
        stat = "binline", alpha = 0.8, draw_baseline = F
    ) +
    geom_point(
        data = toy_nonlin %>%
            group_by(Sex, Treatment) %>%
            summarise(meanEnv = mean(Env), meanResp = mean(Response)),
        aes(x = meanResp, y = meanEnv),
        size = 4
    ) +
    geom_line(
        data = toy_nonlin %>%
            group_by(Sex, Treatment) %>%
            summarise(meanEnv = mean(Env), meanResp = mean(Response)),
        aes(x = meanResp, y = meanEnv, linetype = Treatment),
        linewidth = 1
    ) +
    geom_errorbarh(
        data = toy_nonlin %>%
            group_by(Sex, Treatment) %>%
            summarise(meanEnv = mean(Env), meanResp = mean(Response),
            minResp = mean(Response) - sd(Response),
            maxResp = mean(Response) + sd(Response)),
        aes(y = meanEnv, xmin = minResp, xmax = maxResp),
        linewidth = 1.5, height = 0
    ) +
    # stat_summary(data = toy_nonlin %>%
    #         group_by(Sex, Treatment) %>%
    #         mutate(meanEnv = mean(Env)),
    #     fun = "mean", geom = "line", linewidth = 1,
    #     aes(x = Response, group = meanEnv, y = meanEnv, linetype = Treatment)
    # ) +
    coord_flip() + ylab("Manipulated environmental variable") +
    theme_ridges() +
    theme_bw() + theme(text = element_text(size = 16))

`summarise()` has grouped output by 'Sex'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'Sex'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'Sex'. You can override using the `.groups` argument.
`stat_binline()` using `bins = 30`. Pick better value with `binwidth`.